{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calibration of a Heston / Hull-White model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook is a direct transcription of a QuantLib test case. To test the accuracy of a Heston / Hull-White model calibration, we first generate a data set of option prices, assuming a particular Heston/HW model, then calibrate this model of the option data, and verify that we recover the original model parameters. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set-up"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A  large number of functions is needed to run the calibration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from quantlib.math.hestonhwcorrelationconstraint import (\n",
    "    HestonHullWhiteCorrelationConstraint)\n",
    "\n",
    "from quantlib.time.api import (\n",
    "    Period, Months, March, Date,\n",
    "    Actual365Fixed, TARGET\n",
    ")\n",
    "\n",
    "from quantlib.settings import Settings\n",
    "\n",
    "from quantlib.math.optimization import LevenbergMarquardt, EndCriteria\n",
    "\n",
    "from quantlib.quotes import SimpleQuote\n",
    "from quantlib.termstructures.yields.flat_forward import FlatForward\n",
    "\n",
    "from quantlib.models.api import (\n",
    "    HestonModelHelper, HestonModel, PriceError\n",
    ")\n",
    "\n",
    "from quantlib.processes.heston_process import HestonProcess\n",
    "\n",
    "from quantlib.processes.api import HullWhiteProcess\n",
    "\n",
    "from quantlib.pricingengines.api import (\n",
    "    AnalyticHestonEngine,\n",
    "    FdHestonHullWhiteVanillaEngine)\n",
    "\n",
    "from quantlib.methods.finitedifferences.solvers.fdmbackwardsolver import FdmSchemeDesc\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function creates a flat term structure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def flat_rate(today, forward, daycounter):\n",
    "    return FlatForward(\n",
    "        reference_date=today,\n",
    "        forward=SimpleQuote(forward),\n",
    "        daycounter=daycounter\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The volatility surface is based on a Heston-Hull-White model with\n",
    "the following parameters:\n",
    "    \n",
    "* Hull-White: $a = 0.00883$, $\\sigma = 0.00631$\n",
    "* Heston: $\\nu = 0.12$, $\\kappa = 2.0$, $\\theta = 0.09$, $\\sigma = 0.5$, $\\rho=-0.75$\n",
    "* Equity / short rate correlation: $-0.5$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "dc = Actual365Fixed()\n",
    "calendar = TARGET()\n",
    "todays_date = Date(28, March, 2004)\n",
    "settings = Settings()\n",
    "settings.evaluation_date = todays_date\n",
    "\n",
    "r_ts = flat_rate(todays_date, 0.05, dc)\n",
    "\n",
    "## assuming, that the Hull-White process is already calibrated\n",
    "## on a given set of pure interest rate calibration instruments.\n",
    "\n",
    "hw_process = HullWhiteProcess(r_ts, a=0.00883, sigma=0.00631)\n",
    "\n",
    "q_ts = flat_rate(todays_date, 0.02, dc)\n",
    "s0 = SimpleQuote(100.0)\n",
    "\n",
    "# vol surface\n",
    "\n",
    "strikes    = [50, 75, 90, 100, 110, 125, 150, 200]\n",
    "maturities = [1 / 12., 3 / 12., 0.5, 1.0, 2.0, 3.0, 5.0, 7.5, 10]\n",
    "\n",
    "vol = [\n",
    "0.482627,0.407617,0.366682,0.340110,0.314266,0.280241,0.252471,0.325552,\n",
    "0.464811,0.393336,0.354664,0.329758,0.305668,0.273563,0.244024,0.244886,\n",
    "0.441864,0.375618,0.340464,0.318249,0.297127,0.268839,0.237972,0.225553,\n",
    "0.407506,0.351125,0.322571,0.305173,0.289034,0.267361,0.239315,0.213761,\n",
    "0.366761,0.326166,0.306764,0.295279,0.284765,0.270592,0.250702,0.222928,\n",
    "0.345671,0.314748,0.300259,0.291744,0.283971,0.273475,0.258503,0.235683,\n",
    "0.324512,0.303631,0.293981,0.288338,0.283193,0.276248,0.266271,0.250506,\n",
    "0.311278,0.296340,0.289481,0.285482,0.281840,0.276924,0.269856,0.258609,\n",
    "0.303219,0.291534,0.286187,0.283073,0.280239,0.276414,0.270926,0.262173]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In a first stage, we calibrate a pure Heston model on the data, in order to obtain a good starting point for the Heston/Hull-White calibration, which is much more time consuming."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Heston calibration\n",
      "v0: 0.018030\n",
      "theta: 0.074934\n",
      "kappa: 0.927094\n",
      "sigma: 0.000261\n",
      "rho: -0.134581\n"
     ]
    }
   ],
   "source": [
    "start_v0    = 0.2 * 0.2\n",
    "start_theta = start_v0\n",
    "start_kappa = 0.5\n",
    "start_sigma = 0.25\n",
    "start_rho   = -0.5\n",
    "\n",
    "equityShortRateCorr = -0.5\n",
    "\n",
    "corrConstraint = HestonHullWhiteCorrelationConstraint(equityShortRateCorr)\n",
    "\n",
    "heston_process = HestonProcess(r_ts, q_ts, s0, start_v0, start_kappa,\n",
    "                               start_theta, start_sigma, start_rho)\n",
    "\n",
    "h_model = HestonModel(heston_process)\n",
    "h_engine = AnalyticHestonEngine(h_model)\n",
    "\n",
    "options = []\n",
    "\n",
    "# first calibrate a heston model to get good initial\n",
    "# parameters\n",
    "\n",
    "for i in range(len(maturities)):\n",
    "    maturity = Period(int(maturities[i] * 12.0 + 0.5), Months)\n",
    "    for j, s in enumerate(strikes):\n",
    "\n",
    "        v = SimpleQuote(vol[i * len(strikes) + j])\n",
    "\n",
    "    helper = HestonModelHelper(maturity, calendar, s0.value,\n",
    "                                           s, v, r_ts, q_ts,\n",
    "                                           PriceError)\n",
    "\n",
    "    helper.set_pricing_engine(h_engine)\n",
    "\n",
    "    options.append(helper)\n",
    "\n",
    "om = LevenbergMarquardt(1e-6, 1e-8, 1e-8)\n",
    "\n",
    "# Heston model\n",
    "h_model.calibrate(options, om, EndCriteria(400, 40, 1.0e-8, 1.0e-4, 1.0e-8))\n",
    "\n",
    "print(\"Heston calibration\")\n",
    "print(\"v0: %f\" % h_model.v0)\n",
    "print(\"theta: %f\" % h_model.theta)\n",
    "print(\"kappa: %f\" % h_model.kappa)\n",
    "print(\"sigma: %f\" % h_model.sigma)\n",
    "print(\"rho: %f\" % h_model.rho)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The calibrated parameters are now used as starting point for the full Heston/Hull-White calibration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "'quantlib.pricingengines.vanilla.vanilla.FdHestonHu' object has no attribute 'enableMultipleStrikesCaching'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "Input \u001b[0;32mIn [8]\u001b[0m, in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     12\u001b[0m tGrid \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mmax((\u001b[38;5;241m10.0\u001b[39m, maturities[i] \u001b[38;5;241m*\u001b[39m \u001b[38;5;241m10.0\u001b[39m))\n\u001b[1;32m     13\u001b[0m hhw_engine \u001b[38;5;241m=\u001b[39m FdHestonHullWhiteVanillaEngine(\n\u001b[1;32m     14\u001b[0m     hhw_model, hw_process,\n\u001b[1;32m     15\u001b[0m     equityShortRateCorr,\n\u001b[1;32m     16\u001b[0m     tGrid, \u001b[38;5;241m61\u001b[39m, \u001b[38;5;241m13\u001b[39m, \u001b[38;5;241m9\u001b[39m, \u001b[38;5;241m0\u001b[39m, \u001b[38;5;28;01mTrue\u001b[39;00m, FdmSchemeDesc\u001b[38;5;241m.\u001b[39mHundsdorfer())\n\u001b[0;32m---> 18\u001b[0m \u001b[43mhhw_engine\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43menableMultipleStrikesCaching\u001b[49m(strikes)\n\u001b[1;32m     20\u001b[0m maturity \u001b[38;5;241m=\u001b[39m Period(\u001b[38;5;28mint\u001b[39m(maturities[i] \u001b[38;5;241m*\u001b[39m \u001b[38;5;241m12.0\u001b[39m \u001b[38;5;241m+\u001b[39m \u001b[38;5;241m0.5\u001b[39m), Months)\n\u001b[1;32m     22\u001b[0m \u001b[38;5;66;03m# multiple strikes engine works best if the first option\u001b[39;00m\n\u001b[1;32m     23\u001b[0m \u001b[38;5;66;03m# per maturity has the average strike (because the first\u001b[39;00m\n\u001b[1;32m     24\u001b[0m \u001b[38;5;66;03m# option is priced first during the calibration and\u001b[39;00m\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m     27\u001b[0m \n\u001b[1;32m     28\u001b[0m \u001b[38;5;66;03m# list of strikes by distance from moneyness\u001b[39;00m\n",
      "\u001b[0;31mAttributeError\u001b[0m: 'quantlib.pricingengines.vanilla.vanilla.FdHestonHu' object has no attribute 'enableMultipleStrikesCaching'"
     ]
    }
   ],
   "source": [
    "        h_process_2 = HestonProcess(r_ts, q_ts, s0, h_model.v0,\n",
    "                                    h_model.kappa,\n",
    "                                    h_model.theta,\n",
    "                                    h_model.sigma,\n",
    "                                    h_model.rho)\n",
    "\n",
    "        hhw_model = HestonModel(h_process_2)\n",
    "\n",
    "        options = []\n",
    "        for i in range(len(maturities)):\n",
    "\n",
    "            tGrid = np.max((10.0, maturities[i] * 10.0))\n",
    "            hhw_engine = FdHestonHullWhiteVanillaEngine(\n",
    "                hhw_model, hw_process,\n",
    "                equityShortRateCorr,\n",
    "                tGrid, 61, 13, 9, 0, True, FdmSchemeDesc.Hundsdorfer())\n",
    "\n",
    "            hhw_engine.enableMultipleStrikesCaching(strikes)\n",
    "\n",
    "            maturity = Period(int(maturities[i] * 12.0 + 0.5), Months)\n",
    "\n",
    "            # multiple strikes engine works best if the first option\n",
    "            # per maturity has the average strike (because the first\n",
    "            # option is priced first during the calibration and\n",
    "            # the first pricing is used to calculate the prices\n",
    "            # for all strikes\n",
    "\n",
    "            # list of strikes by distance from moneyness\n",
    "\n",
    "            indx = np.argsort(np.abs(np.array(strikes) - s0.value))\n",
    "\n",
    "            for j, tmp in enumerate(indx):\n",
    "                js = indx[j]\n",
    "                s = strikes[js]\n",
    "                v = SimpleQuote(vol[i * len(strikes) + js])\n",
    "                helper = HestonModelHelper(maturity,\n",
    "                                           calendar, s0.value,\n",
    "                                           strikes[js], v, r_ts, q_ts,\n",
    "                                           PriceError)\n",
    "                helper.set_pricing_engine(hhw_engine)\n",
    "                options.append(helper)\n",
    "\n",
    "        vm = LevenbergMarquardt(1e-6, 1e-2, 1e-2)\n",
    "\n",
    "        hhw_model.calibrate(options, vm,\n",
    "                            EndCriteria(400, 40, 1.0e-8, 1.0e-4, 1.0e-8),\n",
    "                            corrConstraint)\n",
    "\n",
    "        print(\"Heston HW calibration with FD engine\")\n",
    "        print(\"v0: %f\" % hhw_model.v0)\n",
    "        print(\"theta: %f\" % hhw_model.theta)\n",
    "        print(\"kappa: %f\" % hhw_model.kappa)\n",
    "        print(\"sigma: %f\" % hhw_model.sigma)\n",
    "        print(\"rho: %f\" % hhw_model.rho)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "from tabulate import tabulate\n",
    "table = {\"Param\": ['v0', 'theta', 'kappa', 'sigma', 'rho'],\n",
    "         \"Estimated\": [hhw_model.v0, hhw_model.theta, hhw_model.kappa, hhw_model.sigma,\n",
    "                       hhw_model.rho],\n",
    "         \"Exact\": [0.12, 0.09, 2.0, 0.5, -0.75]}\n",
    "\n",
    "print tabulate(table, headers='keys', floatfmt='.3f')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
